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ABSTRACT 


Theoretical  seismograms  for  fundamental  mode  Rayleigh  waves  were 
calculated  for  atmospheric  point  sources  over  oceanic  and  over  continental 
earth  models,  as  recorded  at  an  epicentral  distance  of  10,000  km.  Yields 
were  uniformly  distributed  over  the  range  1  kT  to  10  MT,  for  source  altitudes 
in  the  range  0.3  to  92.0  km.  The  earth  structures  used  were  those  of 
Gutenberg  and  of  Anderson  and  ToksBz.  The  source  models  were  point  mass- 
injection  and  energy-injection  sources  at  altitude,  as  well  as  a  distributed 
pressure  pulse  at  the  surface  of  the  earth. 

It  was  found  that:  (1)  as  far  as  Rayleigh  wave  excitation  is  concerned, 
the  mass-injection  and  energy-injection  sources  are  equivalent;  (2)  for  low 
altitudes  the  Rayleigh  wave  excitation  is  independent  of  source  type,  but 
at  intermediate  altitudes  the  surface  overpressure  source  predicts  greater 
amplitudes  than  the  other  two  source  models;  (3)  for  most  altitudes,  the 
energy  coupling  from  the  atmosphere  into  Rayleigh  waves  is  more  efficient 
for  the  continental  earth  structure  than  for  the  oceanic  structure;  (4) 
Rayleigh  wave  amplitude  is  more  sensitive  to  yield  than  to  burst  height; 

(5)  dependence  of  Rayleigh  wave  amplitude  is  less  than  the  cube  root  relation 
for  low-yield  explosions  at  intermediate  altitudes,  but  greater  for  high- 
yield  explosions  at  near-surface  altitudes;  (6)  spectral  splitting  ratios 
do  not  show  a  systematic  variation  with  yield  and  burst  height. 
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INTRODUCTION 


In  thi 3  paper  we  calculate  Rayleigh  wave  excitation  at  teleseismic 
distances  for  a  variety  of  source  models,  earth  structures,  source  yields, 
and  burst  heights,  and  study  the  variation  of  the  Rayleigh  wave  spectrum  and 
waveform  with  these  parameters. 

We  used  Harkrider's  (1964a)  computer  code  for  calculating  the  far-field 
modal  exditation  of  acoustic-gravity  waves  in  an  atmosphere  consisting  of 
isothermal  layers  overlying  a  rigid  halfspace,  and  a  similar  code 
(Harkrider,  1964b)  for  calculating  Rayleigh  wave  propagation  in  a  layered 
solid/ liquid  half space  underlying  a  vacuum.  These  two  codes  have  been  merged 
in  order  to  model  the  coupling  of  atmospheric  and  seismic  wave  propagation 
and  to  calculate  the  generation  of  Rayleigh  waves  from  a  point  source  in 
the  atmosphere,  and  the  latter  code  was  also  used  to  model  a  pressure  source 
distributed  on  the  solid  free  surface. 

It  was  found  that  the  loading  of  the  atmosphere  on  the  solid  earth 
affected  the  Rayleigh  wave  characteristics  by  only  about  0.01%,  i.e.,  the 
ratio  of  atmospheric  to  solid  earth  density. 

Previous  attempts  to  predict  Rayleigh  wave  excitation  by  atmospheric 
sources  were  either  oversimplified  or  else  heuristic  and  ad  hoc  in  nature. 
Toksoz  and  Ben-Menahem  (1965)  assumed  a  point  source  in  a  homogeneous  liquid 
halfspace  overlying  a  homogeneous  solid  half space;  they  made  use  of  an 
equation  of  Cagniard  (1962)  and  predicted  that  the  effect  of  burst  height 
on  Rayleigh  wave  excitation  should  be  negligible  in  this  situation.  Nickel 
and  Whitaker  (1971)  used  an  ingenious  synthesis  technique  to  approximate  a 
layered  solid  earth  model.  We  discuss  below  the  validity  of  the  assumptions 
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involved  in  their  technique. 

Harkrider  and  Flinn  (1970)  described  the  formulation  for  a  point 
mass-injection  source  within  one  layer  of  the  atmosphere.  In  the  work 
reported  here  we  modified  the  codes  Harkrider  and  Flinn  used  to  include 
Pierce's  (1968)  energy-injection  source  as  well  as  a  distributed  surface 
overpressure  source. 

The  solutions  to  an  expanding  ring  and  an  expanding  disk  surface  load 
on  an  isotropic  solid  half  space  have  been  obtained  by  Gakenheimer  and  Miklowitz 
(1969),  Gakenheimer  (1969),  and  Gakenheimer  (1971)  among  others.  Rayleigh 
waves  and  body  waves  due  to  an  atmosphere  explosion  modeled  by  an  equivalent 
moving  normal  force  applied  to  a  multilayered  half  space  have  been  formu¬ 
lated  by  Hudson  (1969)  and  computed  for  an  equivalent  stationary  point  source 
at  the  free  surface  by  Douglas  et_  al..  (1972) . 

The  codes  used  in  the  present  work  are  described  by  Kerr  (1971) . 

Source  scaling  was  done  by  standard  methods  (Harkrider  and  Flinn,  1970). 

Two  earth  models  were  used  in  the  present  work:  the  Gutenberg  continental 
model  (which  has  a  low-velocity  zone;  see  Ben-Menahem  and  Harkrider,  1964, 
p.  2610)  and  the  Anderson-Toksoz  oceanic  model  (Harkrider  and  Anderson, 

1966,  p.  2970).  Earth  and  atmospheric  layering  characteristics  were  as 
usual  assumed  to  be  invariant  between  source  and  receiver. 

I 

Calculations  were  made  for  a  variety  of  yields  and  burst  heights  for 
both  earth  models  and  for  all  three  source  models.  Waveforms  and  spectra 
were  calculated  fcr  yields  of  1  kT  to  10  MT  and  burst  heights  of  0.3  to 


THEORETICAL  NUMERICAL  METHODS 


Energy  Injection  source  in  the  atmosphere 

The  formulation  for  Rayleigh  waves  generated  by  an  explosion  in  a 
gravitating  atmosphere  coupled  to  the  solid  earth  by  a  mass  injection 
source  was  described  by  Harkrider  and  Flinn  (1970).  We  have  modified  our 


calculations  by  using  the  Pierce  energy  injection  source  (Pierce,  1968 
and  Pierce  et  al.,  1970)  which  was  used  by  them  to  calculate  acoustic- 


gravity  waves  from  an  explosion  in  the  atmosphere. 
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For  an  energy  source  at  (r,  2)  =  (0,  D)  (2  positive  upward)  in  an 

infinite  isothermal  gravitating  ideal  gas,  the  outgoing  time-transformed 

parcel  pressure  <p  >  can  be  written  in  integral  form  as 
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Following  Harkrider  and  Flinn  (1970)  ,  we  find  that 
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Ag  is  the  atmosphere  product  matrix  evaluated  at  the  source  altitude,  and  A 
is  the  same  product  matrix  evaluated  at  the  top  of  the  atmosphere,  k^  is 
the  j'th  root  of  the  Rayleigh  wave  dispersion  equation. 

It  should  be  noted  that  substitution  of  equation  (2)  into  the  source 
relations  of  Harkrider  (1964b)  yields  the  same  acoustic-gravity  wave  over¬ 
pressure  relations  as  those  of  Pierce  et  al.  (1971). 


Rayleigh  wave  displacements  from  sea-level  overpressures 


For  the  excitation  of  Rayleigh  waves  by  an  overpressure  specified  at 
sea  level,  we  model  the  earth  by  a  multilayered  elastic  half  space  with  an 
optional  layered  gravitating  liquid  at  the  surface.  We  formulate  the 
problem  for  an  azimuthally  symmetric  overpressure  source  [pg(r,t)]  applied 
to  the  surface  of  the  earth  model. 

The  Fourier  time- transformed  source 


<pg(r)>  = 


IPg(r,t)]e 


-iwt 


dt 
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can  because  of  symmetry  be  written  as 


<Ps(r)> 


Ps  d0(k*)  dk 


(4) 


where 
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<p  (r)>  J  (kr)r  dr 
s  o 


The  twice  transformed  pressures  and  particle  velocities  at  the  free 
surface  of  the  liquid  layer  for  the  oceanic  model  are  related  to  those  on  the 
ocean  bottom  by 


or 


ps  "  1"  "A21  +  ~  A22 

V  o 

where  the  ocean  layer  matrix  A  can  be  obtained  from  Dorman  (1962)  or  by 
setting  the  gravitational  acceleration  equal  to  zero  in  the  atmospheric  matrices 
of  Press  and  Harkrider  (1962)  and  Harkxider  and  Flinn  (1970). 
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earth  interface  are  related  to  the  elestle  properties  of  the  solid  earth  by 
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-ere  the  q„a„tltles  N,  K,  1,  H.  0,  and  „  are  defined  by  harhrlder  (l964a,  W0) 
Assuming  eo„tlnulty  of  normal  stress  and  vertical  particle  veloclty  at  the 
liquid-solid  interface  we  obtain  by  combining  equations  (5)  and  (6): 
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Tot  the  overpressures  applied  to  the  surface  of  the  continental  model, 
the  result  lg  similar  to  the  oceanic  model: 
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where 
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This  result  could  also  be  obtained  by  letting  the  thickness  of  the  water 
layer  approach  zero  in  equation  (8)  ,  since  k^  -*•  1  and  ^  0  in  thic 
limit. 

Evaluating  the  residue  contribution  of  integrals  (8)  and  (9),  we  obtain 
the  solutions  for  the  vertical  displacement  spectra  associated  with  the  j'th 
mode  Rayleigh  waves.  For  the  oceanic  model  the  solution  is 
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and  for  the  continental  model 


integrate  the  solution  for  a  point  pressure  source  over  the  desired  source 
area.  The  point  pressure  source  it  defined  here  as  pg  =  k/2ir,  in  order 
that  sPg(r)>  given  in  equation  (4)  be  a  delta  function.  Using  equation  (9), 
the  solution  for  the  point  source  is  written 
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The  solution  for  a  distributed  surface  source,  <p  (r,9)>,  is  then  given  by 


the  integral  over  the  surface  S  as 


<w  (r,9)>  = 
o 


"  =  J  <Ps(ro*0o)>  GQL  It0>  dSc 


where  we  define 


G(r|  r  )  = 

~o  2  TT 


LH]  j  (RR)  dk 
FR  ° 


R  =  r  -  r 

1  rv  ^s-O 


and  r  is  the  position  vector  in  the  (r,  9)  plane.  An  example  of  this  technique 
it  given  in  Appendix  A  for  a  disc  source,  which  can  be  verified  by  equation  (4). 


Surface  overpressure  models 


The  distributed  sea  level  overpressures  are  assumed  to  result  from  a 
spherical  wave  incident  at  a  rigid  surface  from  a  point  source  at  a  height  D 
above  the  surface.  The  spherical  wave  front  is  assumed  to  expand  with 

velocity  a.  That  is,  the  overpressure  at  a  distance  r  from  ground 


zero  is 
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<Pe(Ro)>  -  _^s% -  exp  {-i[wR  /a  +  2  ten'1  (m/b  )]} 

<“V>Eo 

where  bs  -  1/T^.  From  equations  (13)  and  (14)  we  see  that 
i-mP  „a 

-  as  s  ■. 

os  2,  ”  2  expl-i2  tan“x(u/b  )] 

0)  +b  s 

s 


i 

(14) 


t 

) 
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Using  the  Sommerfeld  integral,  we  can  express  <p  (r)>  in  the  integral 

s 

form  of  equation  (4)  as 


<Pg  (r)>  ■ 


00 

f  e-lvD 

J.  k  J= 


(kr)  dk 


where 


v2-(f)  -X2 


and  thus 


2cop  a  k  ,  . .  -1 

as  s  exp{-i[vD  +  2  tan  (ui/b  )]} 
0  0  s 

v(w  +b  ) 
s' 


We  obtain  the  ocean  bottom  Rayleigh  wave  displacement  in  the  far  field  by 
substituting  equation  (lb)  into  equation  (10)  and  transforming  back  into  the 


time  domain: 


-  <?> 


i  “I  2  ~y(w)  r 

“paSa0kJ  * 

_ J  _  A 

2  2  i  /  o  ® 

A99v(oj  +  b  )  [a  sin  0  ]  ' 
ll  s  '  e  e 


•  exp{i[(ot-k  r+3ir/4-<J>  ]}  du 
J  s 


where 


i  8  vD  +  2tan-±(co/b  ) 


and 


a  =  i  [GN  -  LH] 

a  (3F  /3k) 

o 

and  where  we  have  used  the  asymptotic  expansion  for  large  arguments  of  the 
Hankel  function,  as  well  as  approximation  for  the  earth's  curvature,  as  in 
equation  (3).  The  expression  for  the  continental  model  can  be  obtained  by 
setting  A22  =  1  and  =  o  in  equation  (16). 


Source  overpressure  scaling 

The  scaling  of  the  outgoing  wave  (equation  13)  for  the  mass  and  energy 
injection  source  models  for  different  yields  and  burst  altitudes  is  done 
using  the  formulas  in  Glasstone  (1962)  Using  values  obtained  from  Glasstone 
(1957)  for  a  1  kT  typical  air-burst  nuclear  explosion,  the  scale  distance 
a  ,  peak  overpressure  p  ,  and  positive  phase  duration  T  are  given  by: 

S  dS  • 3S 


a  =  E.a 
s  *  o 


as 


Ps(D) 


p 

L  o  J 


ao 


T  =  C 

+as 


q 

0 

C  _ 

W 

'Po 

a 

T . 

+ao 

4  - 

P  (D) 

L  s  J 

1/3 


Here,  W  is  the  kT  equivalent  yield,  Pq  and  aQ  are  the  ambient  pressure  and 
acoustic  velocity  at  sea  level,  and  P  and  a  are  at  the  burst  altitude. 
Except  for  the  low  burst  altitudes,  this  scaling  at  the  source  altitude  is 
not  appropriate  for  a  surface  applied  source  where  we  need  the  sea-level  peak 
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overpressure,  unlike  Harkrider  and  Flinn  (1970)  where  the  scaling  of  the 
outgoing  wave  is  done  at  the  same  altitude  as  the  source. 

For  nuclear  explosions  In  the  intermediate  altitude  range,  10  <_  D  <_ 
100  km,  we  use  the  weak  shock  overpressure  relations  derived  by  Murphy 
(1972)  : 


pa  =  0.2  b  [P  (D) 
as  s  s 


,1/2  xe^ 
oJ  b(x) 


(17) 


T 


+as 


.0857  - 
a 


xb(x) 

(1  +  x  +  n) 


b(x)  =  [1-2(1  +  x  +  n)eXE1(x)]1/2 

where  E^(x)  is  the  exponential  integral,  h  is  the  scale  height  of  an 
exponential  atmosphere,  and  x  and  n  are  matching  parameters  for  a  shock 
wave  calculation  in  a  spherically  symmetric  exponential  atmosphere  (Lutzky 
and  Lehto,  1968).  The  matching  parameters  are  given  graphically  and  in 
tabular  form  as  a  function  of  oh  in  Murphy  (1972) .  The  parameter  oh  is  a 
function  of  source  height  and  yield: 


VD) 

29.4  V 


1/3 


where  h  is  in  km,  P  (D)  in  millibars,  and  W  in  kT. 

Equations  (17)  are  the  extrapolation  relations  using  the  Reed-Otterman 
analytic  weak -shock  theory  (Reed,  1959).  As  in  Murphy  (1972),  we  will 
restrict  our  calculations  to  source  heights  and  yields  for  which  the  ratio 
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of  overpressure  to  ambient  pressure  on  the  ground  is  smaller  than  0.1. 
Isothermal  atmosphere  model 

Since  the  result  for  the  atmosphere/solid-earth  system  (equation  3) 
is  similar  to  the  distributed  overpressure  result  (aquation  16),  we  will 
compare  the  simple  model  of  an  isothermal  gravitating  atmosphere  with  the 
sea-levei  distributed  source. 

For  a  single-layer  atmosphere  of  thickness  H, 


ip  e 


-AH 


■  6_ 

3 

01 


sin(vH) 

•v- 


—AH  _ _  _ 

A22  *  e  [cos(vH)  +  f  sin  (vH)] 


where 


f  =  (A  -  g/c2)/v 


and 


6 


4 

U) 


p  is  the  sea  level  ambient  density.  The  A  elements  are  identical  to 
”  8 

the  above,  with  D  replacing  H. 

Thus  in  equation  (3)  for  the  single  layer  system 
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TdT  J  =  **  2AD  ZT  ~  +  if  Ieiy  +  e-1*!) 


vo)  kvH 


e"~[i  +  e-2*  -  if(1  .  e-12vH)J 


where  we  have  used  the  relations  p^D)  .  and  ».  -(H.D)  _ 

Neglecting  the  exponential  terns  involving  the  height  of  the  atmosphere. 


we  have : 


4_  r^i .  fo_ 
■<D>  L  »„  J  ps(D) 


~  3 

VOJ 


exp  [-1.VD  +  i2tan_1f ] 


For  the  range  of  frequencies  anc  phase  velocities  associated  with  Rayleigh 
waves,  the  neglected  exponential  terns  represent  phase  and  group  delays  for 
waves  which  have  travelled  at  least  once  to  the  top  of  the  atnosphere.  As 

H  becomes  larger,  their  contribution  to  the  Rayleigh  wave  arrives  increasing! 

later  than  the  contribution  of  the  direrr  ua„n  r 

direct  wave  from  source  to  ground,  which 

is  represented  in  this  approximation. 

Harbrider  and  Flinn  (1970,  found  that  the  ensphere  had  little  numerics 

effect  on  the  medium  response  Aei  i.e.,  for  the  ocean  system,  A,  and  for 
the  continent  A  =  A 

6  R'  ThUS’  COmparin8  equation  (3)  with  the  continent 
solution  cortesponding  to  equation  (16),  the  integrands,  differ  by  a  factor  S 
for  the  surface  and  Se  for  the  energy  inlection  source,  where  “ 


Sp  =  p  a  e, 
R  as  s  v 


-VD 


and 


S  =  p  a 
e  ao  s  _ 


-vD 


3i2tan-1f 
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after  approximating  6  in  the  Rayleigh  wave  (u,k)  range  by 


6 


4 

tii 


-0) 


4 


For  g  =  o  the  atmosphere  is  no  longer  exponential  and  we  have  f  =  o 
and  v  —  v .  The  two  factors  are  now  equal  since  p  and  p  differ  only  in 
the  ratio  of  ambient  pressures  at  the  source  altitude  and  sea  level.  If  we 
use  the  energy  injection  source  as  our  standard  because  of  its  success  in 
acoustic  gravity  wave  prediction,  the  correct  applied  overpressure,  p  , 

clS 

for  equation  (13)  should  be  the  sea  leve]  value,  p  . 

rao 


NUMERICAL  RESULTS 


Earth  and  atmosphere  models 

Two  earth  models  were  used  throughout  this  study.  The  Gutenberg 
continental  model  was  taken  from  Ben-Menahem  and  Harkrider  (1964,  p.  2610) 
and  the  Anderson-Toksoz  oceanic  model  was  taken  from  Harkrider  and  Anderson 
(1966,  p.  2970).  The  velocities  and  densities  of  these  models  are  shown  in 
Figure  1. 

The  atmosphere  wa3  represented  by  isothermal  layers  comprising  the 
standard  ARDC  atmosphere  (Wares  et  al.,  1960)  and  is  extensively  described 
by  Press  and  Harkrider  (1962)  and  Harkrider  (1964).  The  temperature  profile 
for  the  atmospheric  model  is  shown  in  Figure  2.  The  velocities  c  and  U,  and 
the  medium  responses  A i  for  the  atmosphere-continental  model  and  for  the 
atmosphere— oceanic  model,  as  calculated  by  the  atmosphere— earth  coupled 
system  dispersion  program  HASH,  are  tabulated  in  Table  I. 

Misleading  values  of  Rayleigh  dispersion  for  the  coupled  system  are 
obtained  when  the  Rayleigh  wave  roots 


are  apparently  Intersected  by  the  roots  of  the  atmospheric  acoustic  modes; 

these  occur  along  lines  of  nearly  infinite  slope  in  this  region  of  the 

velocity-period  plane.  Figure  3  shows  the  structure  of  the  phase  velocity 

curves  at  one  such  point  of  intersection,  and  Figure  4  shows  the  group 

velocities  calculated  for  the  roots  shown  in  Figure  3.  The  normalized 

atmospheric  eigen-function  w  /w  is  additional  evidence  for  this  phenomenon: 

mo  r 

Figure  5  shows  w  /w  at  period  41.875  seconds  near  a  point  of  intersection, 
and  T+1  and  T_^  just  0.3125  seconds  above  and  below  T  .  The  periods  chosen 
for  the  Rayleigh  wave  synthesis  described  here  were  required  to  have  appropriate 
kinetic  energy  densities  in  the  bottom  layer  of  the  atmosphere.  This  cri¬ 
terion  resulted  in  the  best  velocity  and  medium  response  values  for  our 
models.  The  group  velocities  for  these  chosen  periods  are  listed  in  Table  I. 

Anelastic  and  scattering  attenuation  coefficient 

The  amplitude  decay  coefficients  -^(to)  used  here  were  taken  from  a  study 
by  Alewine  (1972) .  These  values  differ  appreciably  from  those  used  in 

A 

Harkrider  and  Flinn  (1970),  YB_^(<*0  >  which  were  obtained  from  a  linear 
extrapolation  of  the  long-period  observations  of  Ben-Menahem  (1965). 

A 

The  Y^(“)  values  are  based  on  observations  by  Tryggvason  and  by 
Gutenberg  and  Richter  (Alewine,  1972)  for  continental  paths.  We  have  used 

A 

the  Ya(w)  in  the  Rayleigh  wave  calculations  for  both  continental  and  oceanic 

A 

models.  If  other  coefficients  y^Ou)  are  preferred,  then  the  peak  Rayleigh 
wave  amplitudes  reported  here  could  be  corrected  by 

I 

a  a  2tt 

exp[y  (to  )  -  y  (to  )]  r  where  to  =  — 

A  m  o  m  m  T 

m 
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and  Tm  is  the  period  of  the  peak  amplitude.  If  the  spectral  decay  factor 

with  the  coefficients  Y0(a>)  causes  the  peak  amplitude  Lo  occur  at  a  different 

time  within  the  Rayleigh  wave  train,  corresponding  to  a  different  group 

arrival,  then  the  best  correction  procedure  (the  only  procedure  if  no  wave 

form  is  given)  is  to  redo  the  synthesis  calculation. 

As  an  illustration  of  these  corrections  we  have  calculated  waveforms 

for  a  100  kT  energy  injection  source  at  an  altitude  of  18.6  km  and  epicentral 

distance  of  10,000  km  for  ^(m)  ,  Y>),  and  tf«)  =  0.  The  spectral  decay 

factors  for  a  10,000  km  path  are  shown  in  Figure  6.  We  see  in  Figure  7  that 

the  peak  amplitudes  for  ^  and  ;  occur  at  T  -  21  sec,  but  for  J  at 

m  B-M 

Ti"  “  46  SeC'  Table  11  shou3  a  comparison  between  the  measured  amplitudes  and 

the  amplitudes  oalculated  by  mating  corrections  to  the  measured  y  wave- 

A 

form. 

The  poor  comparison  at  21  sec  for  yq-  0  (123  vs  173)  might  be  the  result 
of  an  anomalous  Airy  phase  amplitude  for  this  synthesis,  since  these  compari¬ 
sons  were  similar  to  those  made  for  other  source  models.  The  peak  amplitudes 

could  also  be  calculated  for  different  epicentral  distances  r'  by  multiplying 
the  amplitude  at  distance  r  by 

"r  sin(r/a  )  11/2 
exp  [y(w  )(r  -  r')]  — - £— 

r'sin(r'/a  ) 
e 

This  correction  is  approximate  and  depends  on  the  amount  of  dispersion 
of  the  wave  train. 
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Results  for  low  source  alti tudes 


The  theoretical  seismograms  presented  by  harkridcr  and  Fllnn  (1970) 
were  recalculated  to  correct  for  a  missing  factor  of  (2tt)_1  and  an 
incorrect  altitude  scaling  for  the  ocean  model.  Examples  of  seismograms  arc 
shown  in  Figure  6  along  with  those  calculated  for  the  energy  injection  source 
(equation  3)  and  the  analytical  surface  overpressure  Glasstone  source 
(equation  16).  For  these  low  source  altitudes,  as  well  as  for  higher  ones, 
the  seismograms  for  the  point  mass  and  energy  sources  are  practically 
indistinguishable.  We  have  chosen  the  energy  source  as  the  standard  because  it 
gives  better  results  for  the  acoustic-gravity  waves.  Table  III  shows  a 
comparison  of  the  peak  amplitudes  for  the  various  low  altitude  source  model.. 

Intermediate  source  altitude  results 

For  source  altitudes  in  the  range  10-1Q00  km  and  yields  of  7.0,  100.  and 
1000  kT,  we  compared  the  Glass  tone-scaled  point  energy  source  (equation  3)  and 
the  Murphy  (1972)  scaled  surface  overpressure  source  (equation  16). 

Typical  waveforms  generated  by  these  source  nodels  are  shown  in  Figure  9. 
The  smooth  spectra  for  the  surface  source  and  the  irregular-looking  spectra 
for  the  source  at  altitude  are  predicted  by  the  source  terms  in  the  integrands 
of  the  respective  equations,  and  demonstrate  the  effect  of  multipathing 
in  the  atmosphere.  The  peak  amplitudes  for  these  source  models,  shown 

in  Figure  10,  are  netrly  the  same  for  the  lowest  source  altitudes,  but  at 
higher  altitudes  the  modified  Sachs  scaling  calculated  by  Murphy  (1972) 
predicts  a  larger  Payleigh  wave  amplitude. 
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Uniform  half  space  models 

In  1971,  Nickel  and  Whitaker  presented  a  numerical  technique  in  which 
the  atmospheric  explosion  was  modeled  by  a  distributed  source  on  the  free 
surface  oi  a  solid  half-space.  The  elastic  constants  of  the  homogeneous 
half-space  were  made  to  vary  with  the  frequency  in  order  to  reproduce 
measured  dispersion  data.  For  a  homogeneous  half-space  the  medium  response 
is  given  by  Harkrider  (1970)  as 


r  * 


ar“  - 


*pvr 


2  2  2  2  9 

(y-1)  +  B  [2V_  -a  -  B  j 
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where 


r  * 
a 


1/2 


and 


Y  -  2 


Assuming  a  Poisson  solid  as  in  Nickel  and  Whitaker  with  the  Lame'  constants 


equal,  then  the  observed  phase  velocity  of  Rayleigh  /aves  V  (u>)  determines 

K 

a  frequency  dependent  shear,  8(u),and  compressions!  velocity,  a (u),  given  by 


VR  -  (2-2/  /3)  S  and 


Thus  the  media  response  for  this  fictitious  half-space  can  be  determined  to 
within  a  factor  of  p,  the  density,  which  is  assumed  to  be  independent  of 
frequency  or 
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Uhcre  VR  Is  In  km/sec  and  T  is  in  seconds.  In  figure  11,  we  show  *R  £„t 
various  earth  srruerures  including  an  oceanic  m^el.  If  „e  assume  aerosol 
density  p  between  2.50  and  3.0  there  is  good  agreement  between  the 
vertically  Inhomogeneous  media  responses  and  the  approximate  half-space  models 
of  Nickel  and  Whitaker  using  our  Rayleigh  dispersion  for  a  continent  and 

ocean.  This  mean.  that  the  technique  of  Nickel  and  Whitaker  is  valid  for 

explosions  over  oceans  as  well  as  mnHn»r,.-D  a 

continents,  a  point  questioned  in  Harkrider 

and  Flinn  (1970). 


yield  and  burst  height  diagnostics 

Analysis  of  theoretical  seismograms  to  obtain  a  set  of  measurements 

capable  of  diagnosing  yield  and  height  of  burst  was  attempted  using  the  energy 

injection  source  model  f„r  both  oceanic  and  continental  layered  structures. 

The  peak  amplitude,  period  of  oscillation  at  the  peak  amplitude,  and  the 

surface  wave  magnitude  were  calculated.  Table  IV  ,  at  an  epicencr.1  distance  of 
10.000  km.  yield.  „f  i  kT  to  10  m>  rad  burst  ,ltltudc8  to  „„  km  ^ 

the  seismic  coupling,  defined  by  the  ratio  of  seismic  energy  to  explosive 
yield,  Es/  w  ,  where 

log  E  -  9.4  +  2.14  H  -  0.054  M  2  (Pomeroy,  1963) 

“  s  s 

and 


Ms  "  log(A/T)  +  1.66  logA  -  0.18 

A  is  pcak-to-peak  amplitude  in  my  anJ  A  is  epicentral  distance  in  degrees. 
Figure  11  show,  the  relation  of  M,  to  yield  and  burst  height  for  the 
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continental  and  oceanic  models.  At  most  altitudes  the  continental  M  is 

s 

greater  than  the  oceanic  Ms  for  a  giv».n  yield;  the  exceptions  are  10  km  and 
51.2  km.  At  a  burst  height  of  80  km,  both  models  show  Mg  proportional  to 
2  l0gV;  and  nCar  0  kra*  Ms  is  Proportional  to  logW.  Figure  13  shows  the  rela¬ 
tion  of  the  maximum  waveform  amplitude  (AJ  to  yield  and  burst  height  for  the 
two  models.  Except  in  the  altitude  range  40-60  km,  the  maximum  amplitudes 
for  the  continental  model  are  greater  than  the  oceanic  model  amplitudes.  If 
a  measured  quantity  could  be  found  with  contours  intersecting  the  contours 
shown  in  Figure  13  at  angles  approaching  90°,  then  we  would  be  able  to  deter¬ 
mine  both  yield  and  burst  height  from  the  Rayleigh  wave.  Without  such  a 
quantity  we  could  determine  the  yield  if  the  burst  height  is  known,  but  knowing 
only  the  yield,  the  burst  height  cannot  be  determined  uniquely. 

If  the  yield  and  burst  height  are  known  for  one  of  a  series  of  explosions 
at  the  same  location,  then  the  yields  for  all  others  detonated  at  approximately 
the  same  altitude  can  be  calculated  by  comparing  the  energy  spectra.  Figure  14 
shows  relative  spectral  energies,  defined  by 


where  N  is  number  of  equally  spaced  spectral  samples  A,  ,  and  the  reference 

X 

spectrum  AR  is  calculated  from  the  theoretical  waveform  for  an  explosion  of 
1  kT  at  the  same  height  and  the  same  epicentral  distance  as  the  event  with 
unknown  yield.  In  terms  of  amplitudes  rather  than  energies,  the  following 
conclusions  can  be  drawn  from  Figure  14 j 

1*  Cube  root  scaling  is  observed  for  the  amplitude  spectra  of 
Rayleigh  waves  from  high-altitude  explosions  over  continental  paths. 
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2.  For  intermediate  altitudes  t-ho  ,  1 /o 

uaes,  the  amplitude  spectra  have  a  W1/2 

dependence  tor  low  yields  and  W1/3  for  high  yields. 

3.  For  near-surface  altitudes  th*  j 

•  the  amplitude  spectra  have  a  W 

dependence  for  high  yields. 

Further  analysis  revealed  that  when  the  amplitudes  differ  from  cube  root 
staling,  the  shorter-period  components  are  closer  to  the  W1'3  dependence 
than  the  longer  period  components.  The  , uantlt,  ER  was  computed  without 

smoothing  the  spectra  A  and  A  ,  hut  the  scatter  In  " 

R’  UL  cne  scatter  in  our  results  might  be 

removed  by  smoothing  at  least 

We  examined  the  amplitude  spectral  minima  and  the  phaae  spectra 

to  see  if  there  existed  a  simple  relationship 
to  the  burst  height,  but  none  could  be  found. 

Spectral  splitting,  l.e..  the  ratio  of  energy  in  the  longer-period  end 
of  the  Rayleigh  wave  spectrum  to  energy  In  the  shorter-period  end,  has  been 
suggested  as  a  discriminant  between  under  -round  explosions  and  earthquakes 
(von  Seggem  and  Laatert,  1,70).  The  physical  rationale  is  that  the  two 
source  types  should  be  sufficiently  different  in  nature  that  the  shape  of  the 
radiated  elastic  wave  spectrum  should  be  measurably  different.  Intuitively 
it  seems  reasonable  that  variations  in  yield  and  burst  height  of  atmospheric 
explosions  might  produce  similar  effects,  so  we  studied  the  spectral  distri¬ 
bution  of  energy  in  the  theoretical  Rayleigh  wave  calculations. 

We  define  the  power  spectral  splitting  ratio  by 


p(t2,t) 


T  Pd.T^  ’  V  T  1  t2 


Where  P(a,b) 


i=a  Ai  Snd  the  Ai  are  the  F°urier  amplitude,  coefficients. 
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von  Seggern  and  Lambart  (1970)  showed  that  tha  valuaa  of  obtainad  by 
splitting  tha  apattra  at  I  -  22  aaaonds  aaparatas  tha  population  of  atmospheric 
explosion  sals.ogtams  from  earthquakes  and  undatgtound  axploslon  seismograms 
hatausa  tha  former  contain  lasa  long-period  energy  than  tha  latter.  Putting 

T1  "  10  seconds  and  t2  ’  62  seconds,  we  obtain  the  values  of  <R>22  shown 
In  Figure  15  for  tha  continental  nodal.  The  only  general  relation  we  can 
conclude  from  these  data  is  that  for  tha  higher  burst  altitudes  there  is  a 
yield  threshold  above  which  <R>22  is  constant.  This  relation  is  also  observed 
in  Figure  16,  which  gives  <R>T  for  T  -  15,  20,  25,  30,  and  40  seconds  as  a 
function  of  burst  altitude  for  yields  of  1,  10,  and  10,000  kT.  In  another 

attempt  to  use  <R>T  for  yield  and  burst  height  diagnostics,  we  calculated 
fhe  function 


R(T)  =  log  <r>t] 


and  we  found  that  the  shapes  of  R(T)  were  nearly  independent  of  the  yield 
but  were  distinguishable  for  different  burst  heights.  (Figure  17).  Unfortu¬ 
nately  there  does  not  appear  toJ>e  any  systematic  variation  that  would  permit 
a  burst  height  prediction  from  R(T)  when  the  source  parameters  are  unknown. 


CONCLUSIONS 


Our  theoretical  treatment  of  Rayleigh  waves  generated  by  atmospheric 
nuclear  explosions  contains  a  number  of  approximations  with  regard  to  source 
spectra  and  scaling,  the  atmospheric  model,  and  the  earth  model.  Improvements 
in  these  assumptions  will  require  comparisons  with  observations  from  sources 
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whose  burst  height  and  yield  are  known.  The  choice  of  peak  overpressure 
scaling  depends  on  such  considerations. 

Having  analyzed  our  theoretical  data  base  to  find  source  diagnostics, 
we  found  that  source  yield  may  be  estimated  if  the  source  altitude  is 
approximately  known.  We  expect  the  variation  of  yield  to  result  in  a 
continuous  change  of  the  waveforms  and  spectra.  In  the  case  of  source 
altitude  we  are  dealing  with  higher  modes  of  the  atmosphere  system  due  to 
the  large  Rayleigh  phase  velocities  compared  to  the  acoustic  velocities 
of  the  atmosphere.  This  means  there  are  many  nodal  points  in  the  atmosphere 
eigenfunction  with  respect  to  altitude  and  frequency.  Since  this  eigen¬ 
function  is  a  factor  in  the  source  altitude  term  [see  equation  (3)]  the 
variations  in  the  altitude  resulted  in  nearly  discontinuous  changes  in 
spectra.  We  conclude  that  the  yield  diagnostics  may  not  be  useful  for 
determining  source  altitude. 

The  seismic  coupling  is  generally  greater  for  explosions  over  continents 
than  over  oceans;  the  difference  may  be  as  much  as  0.5  magnitude  units.  The 
major  differences  resulting  from  using  different  earth  models  are  seen  in 
the  waveforms,  which  are  governed  largely  by  the  group  velocity  dispersion. 
Since  most  travel  paths  are  a  combination  of  oceanic  and  continental  segments, 
the  effective  dispersion  will  be  a  weighted  average  of  the  velocities  used 
in  computing  these  waveforms.  The  real  earth  paths  will  also  cause  changes 
in  the  amplitude  spectra  due  to  mode  conversions  and  reflections  where  the 
structure  is  changing.  In  consideration  of  these  effects  (which  we  have  not 
included  in  our  analysis)  we  conclude  that  the  best  estimates  of  source 
altitude  and  yield  will  be  made  by  comparing  seismograms  from  explosions 
having  common  epicenters  and  travel  paths  such  as  the  technique  used  by 
Toksoz  and  Ben-Menahem  (1964) . 
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Figure  6.  Spectral  decay  factors  cxp(-Y,\A)  and  exp  (-yb-Ma)  for 
A  =  10000  km  cpiccntral  distance. 
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Figure  7.  Waveforms  for  different  spectral  decay  factors  described 
in  text.  Source  altitude  is  18.6  km,  yield  is  100  kT,  and  cpiccntral 
distance  is  10000  km. 
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Figure  8a.  Seismograms  for  1  kT  yield  and  A  =  10000  km  for 
different  source  types,  source  altitudes,  and  earth  models. 
Continental  model,  source  at  0.305  km.  In  each  case  the  upper 
seismogram  for  a  mass-injection  source,  the  middle  seismogram 
for  an  energy-injection  source,  and  the  lower  seismogram  for 
a  surface  overpressure  source. 
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Figure  8c.  Seismograms  for  1  kT  yield  and  A  =  10000  km  for 
different  source  types,  source  altitudes,  and  earth  models 
Oceanic  model,  source  at  0.305  km.  In  each  case  the  upneJ 
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gram  or  a  surface  overpressure  source. 


I 


different °source°types  1  kT  yield  and  A  =  10000  km  for 
Oceanic  models,  soulrce’at  T.83  km^fn^ \3nd  earth  models 
seismogram  for  a  mass-ini ectinn  ’  In  each  case  the  upper  * 
gram  for  an  energy-iniecti on  source,  the  middle  seismo- 
8™  f°r  =  surface  ovn^““’°Usr^cea?d  ^  lmrer  «*»- 


41 


5  0  km/sec 


S  O  km/ssc 


.25  m|i 


100  sec  2.2  km/sec 


T  100 


i  10 


10.00  20.00  3.00  40.00  50.00  60.00 

PERIOD  (sec) 

Figure  9b.  Seismograms  and  amplitude  spectra  for  intermediate 
source  altitudes.  Continental  earth  model,  for  epicentral 
distance  A  =  10000  km.  The  upper  figures  show  results  for  an 
energy  injection  source  and  the  lower  fi^  ires  for  a  surface 
overpressure  source,  100  kT  at  18.6  km  altitude. 
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Figure  12.  Ms  as  a  function  of  source  altitude  with  yield  as  a 
parameter,  for  continental  and  oceanic  structures. 


ure  13.  Maximum  Rayleigh  wave  amplitude  contoured  as  a  functi 
yield  and  burst  height,  for  continental  and  oceanic  structure 
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igure  16 .  Spectral  power  ratio  <R>22  as  a  function  of  source 
ltitude ,  for  T  =  15  to  40  seconds,  and  for  three  source  yields, 
ontinental  earth  model. 
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TABLE  2 

Peak  amplitude  corrections  for 
different  attenuation  factors 


Attenuation 

coefficient 

T  (sec) 
m 

Synthesized 
Waveform  Peak 
Amplitude 

Amplitudes 

calculated  from 

Y.  waveform 
A 

o 

11 

o 

21 

173 

123 

46 

74 

75 

Y 

A 

21 

64 

46 

39 

^B-M 

21 

5 

5 

46 

17 

18 
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TABLE  3 


Peak  Rayleigh  wave  amplitudes  (mp)  for 
1  kT  yield  at  low  altitudes 

Continental  Structure  Oceanic  Structure 


Source  height 
(km) 

Mass 

source 

Energy 

source 

Surface 

Overpressure 

source 

Mass 

source 

Ener  gy 
source 

Surface 

Overpressure 

source 

0.31 

0.51 

0.49 

1.06 

0.49 

0.44 

0.71 

1.83 

1.06 

1.08 

1.09 

0.58 

0.62 

0.69 

3.66 

1.04 

0.97 

0.99 

0.88 

0.86 

0.71 

4.88 

1.03 

1.08 

1.05 

0.45 

0.52 

0.72 

y 
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TABLE  IV  Cont’d 


I 

I 

I 

I 


28.1  4.28  .000925 


Continental  Model 


APPENDIX  A 


Disc  pressure  source 


Our  assumed  pressure  source  is  constant  over  a  disc  of  radius  a 

and  with  a  transformed  time  history  B(u>),  i.e.. 


B(w), 


<Ps(r,6)>  = 


(A-l) 


Using  the  area  integration  technique,  we  obtain 


2  it  ,a 


<w  (r,0)>  = 


/  l 


Gf£l£o>  rodrod0l 


( 

■5,w  / 


.  2tt  -a 


!  /  - 


o(VR)ro  drQ  d6o 


where 


p.^  =  +  r^  -  2rr  cos  (6-0  1 


From  the  addition  theorem  of  Bessel  functions, 


J  (kR)  =  T  e  J  (kr  )  J  (kr)  cos(m(6-0  )) 
r,  '  m  m  o  m  o 


where  the  Neuman  factor  e  is  defined  as 

m 


1 ,  m  =  o 


2  ,  m  ^  o 


Since 


cos[m(6-eo)]  =  cos(m9)cos(m0o)  -  sin(m6)sin(meo) 


the  integration  over  0  yields 


<w  (r ,0) >  =  -iB(ui) 


f  [GN  -  LHl  ....  / 

0  -  Jo  <kr>  J  Jo  <kr  )rQ  drQ  dk 


where  we  have  used  the  relations 


cos(m0)  d0  = 


o  ,  m 


2tj  ,  m  -  o 


sin(m@)  d@  =  o  for  all  m. 


/ 


Since 


t  Ji  <k*> 


our  solution  for  the  disc  source  is: 


<w  (r,6)>  =  -iB(u)a 
o 


[GN  -  LH 


1 

J 


J^(ka) 

k 


J  (kr)  dk 
o 


(A-2) 


This  solution  can  be  verified  using  equations  (4)  and  (9) .  The  disc  source 
can  be  expressed  in  integral  form  as 


<p  (r,0)>  =  B(u)a 
s 


J,  (ka)  J  (kr)  dk 
1  o 


and  comparison  with  equation  (4)  implies  that 

P  =  B(oj)  a  Jt  (ka) 
s  i 


which  substituted  in  equation  (9)  gives  the  same  result  as  equation  (A-2). 
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